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Abstract 

The need for efficient and controlled expansion of cell populations is paramount in tissue engineering. Hollow fibre 
bioreactors (HFBs) have the potential to meet this need, but only with improved understanding of how operating 
conditions and cell seeding strategy affect cell proliferation in the bioreactor. This study is designed to assess the effects of 
two key operating parameters (the flow rate of culture medium into the fibre lumen and the fluid pressure imposed at the 
lumen outlet), together with the cell seeding distribution, on cell population growth in a single-fibre HFB. This is achieved 
using mathematical modelling and numerical methods to simulate the growth of cell aggregates along the outer surface of 
the fibre in response to the local oxygen concentration and fluid shear stress. The oxygen delivery to the cell aggregates 
and the fluid shear stress increase as the flow rate and pressure imposed at the lumen outlet are increased. Although the 
increased oxygen delivery promotes growth, the higher fluid shear stress can lead to cell death. For a given cell type and 
initial aggregate distribution, the operating parameters that give the most rapid overall growth can be identified from 
simulations. For example, when aggregates of rat cardiomyocytes that can tolerate shear stresses of up to 0.05 Pa are evenly 
distributed along the fibre, the inlet flow rate and outlet pressure that maximise the overall growth rate are predicted to be 
in the ranges 2.75 x 10^'m-s^' to 3 x lO^^m^s^' (equivalent to 2.07mlmin^' to 2.26mlmin^') and 1.077 x 10' Pa to 
1.083 x 10' Pa (or 15.6 psi to 15.7 psi) respectively. The combined effects of the seeding distribution and flow on the growth 
are also investigated and the optimal conditions for growth found to depend on the shear tolerance and oxygen demands 
of the cells. 



Citation: Chapman LAC, Shipley RJ, Whiteley JP, Ellis MJ, Byrne HM, et al. (2014) Optimising Cell Aggregate Expansion in a Perfused Hollow Fibre Bioreactor via 
Mathematical Modelling. PLoS ONE 9(8): el05813. doi:10.1371/journal.pone.0105813 

Editor: Ben D. MacArthur, University of Southampton, United Kingdom 

Received May 26, 2014; Accepted July 24, 2014; Published August 26, 2014 

Copyright: © 2014 Chapman et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: The work is primarily supported by an EPSRC Life Sciences Interface Doctoral Training Centre grant from the University of Oxford {EP/F500394/1 ). HMB 
is supported by King Abdullah University of Science and Technology, Saudi Arabia (Award No. KUK-Cl -01 3-04). SLW is funded by the EPSRC through an Advanced 
Research Fellowship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 

Competing interests: The authors have declared that no competing interests exist. 

* Email: waters@maths.ox.ac.uk 



Introduction 

In vitro tissue engineering that uses cells extracted from the 
patient to grow tissues in the laboratory offers a viable alternative 
to donor transplants for replacing tissue that has been damaged or 
lost. However, to grow substitutes of clinically relevant dimensions 
it is necessary first to expand the population of cells from the 
patient. Hollow fibre bioreactors (HFBs) are promising cell 
expansion devices: they consist of a glass module containing a 
single or multiple hoUow porous biodegradable polymer fibre(s) 
onto (or around) which the cells are seeded. In this study, we 
consider the single-fibre HFB shown in Figure lA. The outer 
surface of the fibre provides a large area (relative to the volume of 
the device) for cell proliferation. The flow rate of culture medium 
into the lumen, g,„, is prescribed and the relative fluid fluxes 
through the fibre wall (or membrane) and out of the lumen are 
controlled by a back pressure, Pi^out, applied at the lumen outlet. 
This enables the nutrient and fluid shear stress environment of the 
cell population to be controlled to encourage functional growth. 
Two outiets from the extra-capUlary space (ECS), the space 
around the fibres, can be closed or opened to promote flow 



through the membrane and improve nutrient delivery to the cells 
and removal of waste products [1]. 

Despite these advantages, the combined effects of the flow 
conditions, nutrient transport and initial cell seeding distribution 
on cell population expansion are not well understood. Increasing 
the flow through the membrane improves nutrient delivery to the 
cells but also increases the fluid shear stress they experience. Whilst 
higher nutrient concentration promotes cell proliferation and 
higher shear stress can stimulate proliferation, excess shear stress 
can lead to cell death. Hence it is necessary to find the flow 
parameters and seeding distribution that give the optimal balance 
of nutrient delivery and shear stress for the given cell type. 

Cells are typically seeded onto the surface of the fibre by 
injecting an inoculum via the ECS outlets and fixing the module to 
a drum rotating at 6.6 rpm for 6hrs [2]. This reduces cell damage 
and enhances cell adhesion compared to pumping the inoculum 
through the ECS. However, the cell distribution over the surface 
after culturing is often uneven (Figure IB) and the yield poor 
compared with that achieved on flat tissue culture plastic [2]. 
During seeding, cells attach to the membrane individually or in 
small groups, but many cell populations grow predominantly as 
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Figure 1. Hollow fibre bioreactor and cell distribution following culture. A, Photograph of single-fibre HFB module. ECS = extra-capillary 
space. Adapted from [1]. B, Uneven distribution of pZIP osteogenic cells on a poly(lactic-co-glycolic) acid hollow fibre 7 days after seeding. Dark 
patches are methylene-blue-stained cells. Top section is from the outlet end of the fibre, middle section from the middle of the fibre and bottom 
section from the inlet end with the flow going from right to left. Photo reproduced from [2] with permission. 
doi:1 0.1 371/journal.pone.01 0581 3.g001 



aggregates [2-5] . Experimental and modelling studies have sought 
to explain the interactions between the flow of culture medium, 
initial cell distribution and cell population growth in specific 
perfusion bioreactor systems (see [6] for a review). Although 
general relationships have been identified, they typically pertain to 
cells seeded in 3D scaffolds, so have limited applicability to the 
HFB. Indeed, no study has considered the effect of cell seeding on 
cell proliferation in a HFB. While it is possible to mc'asure the 
nutrient concentration globally in HFB experiments, for example 
at the inlet and outlets, there is a lack of detailed spatial 
information on local nutrient levels, fluid shear stress and cell 
distribution inside the bioreactor. 

These issues motivate the development of the following 
mathematical model, in which we consider aggregates of cells 
that an; a single c:ell thick growing along the outer surface of 
the fibre in a HFB. The growth rate of the aggregates depends on 
the local oxygen concentration (as oxygen is considered to be the 
growth-rate-limiting nutrient [1]) and the shear stress exerted on 
them by the flow (following [7,8]). We determine optimal 
operating conditions (lumen inlet flow rate and outlet pressure) 
for growth of the aggregates to confluence for different cell types, 
and investigate the impact of the seeding distribution on the 
aggregate growth rate. 

Numerous experimental and theoretical studies of fluid and 
nutrient transport in HFBs have been carried out (see [9]). 
Theoretical studies have typically focussed on fluid and mass 
transport in a single-fibre unit (referred to as a Krogh cylinder). In 
the lumen these are described by the Stokes and advection- 
diflusion equations respectively, while in the membrane and ECS 
com'ective effects are typically ignored and mass transport is 
assumed to be diffusion-dominated. This is representative of the 
setup with the ECS outiets closed. If present, cells are usually taken 
to be seeded in gel throughout the ECS so that nutrient uptake can 
be modelled by a reaction term in the ECS mass transport 
equation. This reaction term has followed zeroth-order, first-order 
or fuU nonlinear Michaelis-Menten kinetics depending on the 
nutrient demands of the cells, and the nutrient transport equations 
have been solved via appropriate analytical or numerical methods 
(see [9] for a review). Shipley and Waters [1] adopted an analytical 
approach to model fluid and nutrient transport in a single-fibre 
HFB, with cells seeded in gel in the ECS. Their work represents a 
departure from previous studies by considering nutrient advection 



in the membrane and ECS (representative of open ECS outlets 
and/or an applied downstream lumen pressure), using Darcy's 
Law and reaction-advection-diffusion equations to describe the 
fluid flow and mass transport respectively. Growth of the cell 
population and its effect on the permeability of the ECS to fluid 
and nutrient were neglected. This is typical of theoretical studies of 
HFBs, which focus on timescales associated with transport 
processes rather than those associated with cell proliferation, and 
model the ECS as devoid of cells or packed with cells already at 
physiological densities, so that proliferation may be neglected [10]. 
The models of Brotherton and Chau [10] and Mohebbi-Kalhori et 
al. [1 1], in which Monod kinetics were used to describe nutrient- 
dependent cell proliferation in the ECS, are notable exceptions. 
Brotherton and Chau also included feedback of the cell 
proliferation on tlu- fluid flow and nutrient uptake by making 
the ECS fluid permeabihty and nutrient uptake functions of the 
cell density. 

None of these studies, however, have considered the influence of 
fluid shear stress on cell proliferation. Experimental studies have 
investigated the effect of fluid shear stress on the differentiation 
and proliferation of bone cell progenitors (see [12,13] for 
comprehensive reviews), which are commonly ctiltured in HFBs. 
Computational approaches have been used to assess the role of 
shear stress in tissue construct development [8,14-16]. Shakeel et 
al. [8] modelled shear-stress-dependent cell proliferation and 
nutrient transport in a perfusion bioreactor. They adapted a 
functional form proposed by O'Dea et al. [7] to describe elevated 
proliferation at intermediate shear stress levels, and used a similar 
function to describe faster nutrient uptake at intermediate shear 
stresses. Korin et al. [15] combined modelling and expc'rim(;nts to 
investigate the effect of oxygen transport and shear stress on 
human foreskin fibroblast proliferation in a microchannel biore- 
actor. Cell proliferation was governed by a logistic growth law and 
oxygen uptake at the cell surface described via Michaelis-Menten 
kinetics, which were independent of the changing cell density. 

The model presented here differs from previous studies of HFBs 
in a number of key respects. We assume that cells are attached to 
the outer surface of the fibre rather than seeded throughout the 
ECS, so the model has greater relevance to cell population 
expansion. Cell proliferation and death are modelled phenome- 
nologicaUy by nutrient- and shear-stress-dependent elongation and 
shortening of the cell aggregates. This enables us to determine the 
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impact of changes in tlie operating parameters and seeding 
conditions on tlie growth rate of tlie aggregates, rather than just on 
the fluid and mass transport through the bioreactor. Additionally, 
we account for feedback on the flow and nutrient distribution from 
cell proliferation and death (via changes in the permeability of the 
outer membrane surface and oxygen uptake with aggregate 
growth). 

The structure of the paper is as follows. In the Methods section, 
we present the governing equations for fluid and oxygen transport 
and cell aggregate growth, and identify values for the transport 
parameters. We briefly describe how the model is parameterised 
for specific cell types, simplified and solved numerically (further 
information is given in File SI). In the Results section we perform 
a parameter sensitivity analysis, assessing the impact of changing 
the inlet flow rate, outlet pressure and initial aggregate distribution 
on the aggregate growth, and identify optimal operating condi- 
tions that minimise time to confluence for a specific cell type. We 
summarise our findings and explain how they may be used to 
maximise the cell yield for different cell types in the Discussion 
section. 

Methods 

Model Setup 

A schematic of the HFB is shown in Figure 2 A. We consider a 
representative 2D cross-section as shown in Figure 2B. We assume 
that the ECS outlet nearest the lumen inlet is closed while that 
nearest the lumen outiet is open. Culture medium is pumped into 
the lumen at a prescribed rate, Qm, and flows out of the bioreactor 
via the lumen or the open ECS outiet. The pressure imposed at the 
lumen outiet, P/^out, determines the ratio of flow through the 
lumen to that through the membrane into the ECS. All external 
boundaries other than the inlet and outlets are impermeable to 
fluid and oxygen. The oxygen concentration in the culture 
medium pumped into the inlet, C,„, is constant. Oxygen is 
transported through the bioreactor by advection and diffusion. For 
simplicity, we treat the flow and oxygen distribution as being 
symmetric about the lumen centreline. Cell aggregates are 
distributed along the outer surface of the membrane, and absorb 
oxygen and grow or shrink at a rate that depends on the local 
oxygen concentration and fluid shear stress. This change in 
aggregate length affects the flow and oxygen profiles along the 
membrane through increased/decreased oxygen uptake and 



reduced/increased permeability of the membrane outer surface 
to fluid. 

Governing Equations 

Fluid Transport. We let x denote the axial distance down 
the lumen centreline, with the inlet at x = 0 and the outlet aXx = L, 
and y the height above the lumen centreline. The corresponding 
unit vectors in the x- and j-directions are denoted by e.v and e^. In 
the lumen and ECS, denoted by subscripts / and e respectively, we 
treat the flow as quasi-steady on the timescale of aggregate growth, 
which is much longer than the timescale for fluid transport (days or 
weeks compared to seconds or minutes). We neglect inertial forces 
in the fluid since they are much smaller than viscous forces for 
experimentally-relevant parameter values (this is verified a poster- 
iori — see the Reduced Model section). Hence, we describe the 
fluid flow by the steady incompressible Stokes flow equations 

V-u,=0, ^iV\ = Vpi, i = l,e, (1) 

where U, = (m/,v,), pi and p. are the fluid velocity, pressure and 
dynamic viscosity respectively. In the porous membrane, denoted 
by subscript m, fluid transport is modelled using the incompress- 
ible Darcy flow equations 

k 

V-u„, = 0, u„,= ^Vp„,, (2) 

where Um and Pm are the fluid velocity and pressure (averaged over 
the pore space), and k„, is the membrane permeability. 

At the interface between the lumen and the membrane, y = H, 
we impose continuity of normal velocity and normal stress. We 
make the modelling assumption that all of the lumen fluid stress is 
taken up by the fluid in the membrane as in previous studies [1,17] 
(other constitutive relations could be considered, see [18] for a fuU 
discussion). Hence the continuity conditions at the interface are 

u/-e3, = 0u„,-ej,, e,.-o-/-eJ = e^/(^(7,„-eJ ony = H, (3) 

where (/> is the porosity of the membrane (assumed constant), and 
(7,- (i = l,m,e) are the fluid stress tensors given by 



ECS outlet y 

[/ 




Figure 2. Schematics of hollow fibre bioreactor and model setup. A, Schematic of the laboratory single-fibre HFB module. B, Schematic of 
model setup. Cell aggregates on membrane shown in black. Arrows show direction of fluid flow. Red rectangle in A shows the region that B 
corresponds to. H= lumen radius, //,„= membrane radius, //, = ECS radius, L= length of model domain, 2,„ = lumen inlet flow rate, 
lumen outlet pressure, and P„„„ = atmospheric pressure. 
doi:1 0.1 371/journal.pone.01 0581 3.g002 
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<ri= -pil+K^Ui+'VuJ), i=l,e, am=-PmI- (4) 



Ueex = 0 onx = 0, Hm<y<He. 



(9) 



The cell aggregates on the outer surface of the membrane are 
taken to be infinitesimally thin (sinc:e a topical cell diameter 
(10|im) is small compared with the depth of the membrane 
(200 |im) and the aggregates are one-cell thick). At the outer 
surface of the membrane, 7 = //m, we impose continuity of normal 
fluid velocity but assume there is a jump in the normal stress 
proportional to the fluid flux, so that 

0u„-ev=u<,-ev, 



K(x,t) 



Ue Cy on y=Hm 



(5) 



where K{x,t) represents the permeability of the outer surface of 
the membrane. We assume that the aggregates provide resistance 
to the flow such that K{x,t) has a low value in regions covered by 
cell aggregates. In regions not covered by aggregates, we assume 
that K(x,t) has a higher but finite value, Ki,i, to model a reduction 
in surface permeability caused by extra-ceUular matrix and protein 
deposition by the cells. The discontinuity in K(x,t) at the ends of 
the aggregates is smoothed in our numerical simulations as 
described in Section C.2 of File SI. We suppose that the cell 
density within an aggregate remains constant as the aggregate 
grows/shrinks, so that Kig and Ki,i are fixed. The normal stress 
boundary condition in (5) corresponds to Starling's equation for 
the fluid flux across a thin membrane, i.e. the fluid flux being 
proportional to the pressure difierence across the membrane (see 
Section A.l in File SI). 

At the interfaces between the membrane and the free fluid we 
impose no-slip conditions 

ure^ = <l>u,„e^ ony = H, 



pu„ex = ne-ex on y = Hm, 



(6) 



since Shipley et al. [19] fitted a mathematical model that 
incorporated slip to experimental measurements of fluid distribu- 
tion and found that slip was negligible for these membranes. We 
assume the flow is symmetric about the lumen centreline, so 



u/-e^=0, Vu/-ej,=0 0113^ = 0. 



(7) 



The outer boundary of the ECS is impermeable, so we impose 
no-slip and no-flux conditions there 



Ue = 0 on y = Hg. 



(8) 



Having specified the transverse boundary conditions for the 
fluid flow we now specify the appropriate axial boundary 
conditions. There is no fluid flux out of the ends of the membrane 
or at the upstream end of the ECS 

Umex = 0 onx=0,I,, H<y<H,„, 



The flow rate at the inlet and the normal stress at the lumen 
outlet are prescribed so that 



U/-exlx=odj=2;«. (10) 



ex-oi-^l=-Pi,out onx=L, Q<y<H, (11) 

where the values of 2/„ and Pioui are taken from experiments. 

Although solving full Stokes flow would rerjuire two pointwise 
boundary conditions at the inlet, we only require (10) since we 
exploit the small aspect ratio of th(- fibri; to simplify the flow 
equations (see the Reduced Model section). At the ECS outlet we 
impose continuity of normal stress 



onx = L, H„<y<He, 



(12) 



where Pam is atmospheric pressure (« 1.013 X 10^ Pa or 
14.69 psi). 

Oxygen Transport. We assume that oxygen transport is 
quasi-steady on the timescale of aggregate growth, and governed 
by a steady advection-difiusion equation in each region 



u,-Vc, =Z),V^c,-, i = l,m,e. 



(13) 



where c,- (i = l,e) and Z),- {i = l,e) wet the oxygen concentration (per 
unit volume of fluid) and difiFusivities in the relevant domains, and 
Cm and Dm are the average oxygen concentration and diffusivity in 
the membrane. We assume that the diffiasivities are constant. 

At the lumen-membrane interface, we impose continuity of the 
oxygen concentration and flux 



ci = Cm, Di\ciey = Dm^Cmey on y = H. 



(14) 



At the interface between the membrane and the ECS, we again 
impose continuity of the concentration, but also impose a jump in 
the flux across the aggregates due to ceUular uptake of oxygen 



Cffi Cg, 
R(ce) in an aggregate, 



(DeVce-DmVc„yey = i (15) 
0 outside an aggregate, 

on 3; = H^, 

where the oxygen uptake flux R(Ce) is an increasing saturating 
function of Ce 



Ci/2 + Ce 



with the positive constant F denoting the maximal oxygen uptake 
flux and C1/2 the concentration at which the flux is half-maximal. 

Since the oxygen distribution is taken to be symmetric about the 
lumen centreline, we impose no difiusive flux through the lumen 
centreline 
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Vc;e^ = 0 on 3^ = 0. 



(16) 



The upstream end and outer wall of the ECS and the ends of 
the membrane are impermeable to nutrients, so that 

Vce-e^ = 0 onx=0, H„<y<He, (17) 

Vc„e;c = 0 onx = 0,l, H<y<H„, (18) 

Vce-ey = 0 ony = He. (19) 

The inlet oxygen concentration is fixed, so 

ci = Ciri onx = 0, 0<j<i/. (20) 



Following Shipley et al. [1], we assume that the culture medium 
leaving the lumen and ECS outlets is well-mixed and impose zero- 
diffusive-flux boundary conditions at the outiets, 

Vc; e;c = 0 onx=L, 0<y<H, 



Vceex = 0 onx=L, Hm<y<He. 



(21) 



Cell Aggregate Growth. Net increases/ decreases in the cell 
population are modelled by elongation/shortening of the aggre- 
gates, whose initial distribution and lengths along the membrane 
are prescribed. It is assumed that cell proliferation occurs over the 
whole length of the aggregate, and that when cells in the middle of 
an aggregate divide the neighboring cells move outwards to 
accomodate the new cells, maintaining the cell densit)' as the 
aggregate grows. This is based on evidence from cell migration 
and proliferation assays on 2D substrates that cells distant from a 
growing edge continue to proliferate and exert a mitotic pressure 
on their neighbouring cells [20-23]. When two aggregates come 
into contact, they are assumed to coalesce and are thereafter 
treated as one larger aggregate. If one end of an aggregate reaches 
either end of the membrane (x = 0 or x = L), any further growth is 
taken to occur at its opposite end. It is assumed that when cells die 
they detach from the membrane and are removed with the flow, 
causing the aggregates to shrink. For simplicit)', we assume that the 
eflects of the oxygen concentration and fluid shear stress on the 
cell proliferation rate Gp and death rate Gj, and hence on the net 
rate of aggregate elongation, are multiplicative (see Equation (22)). 
This means that the aggregates cannot grow in either hypoxic or 
excess shear conditions, which would be possible if the effects were 
assumed to be additive (e.g. if the shear stress was high enough to 
kill cells, but the oxygen concentration was also ver^' high and 
promoted proliferation suflficiendy for the net effect to be positive). 
Nevertheless, the modular nature of the model means that 
alternative growth laws (based on different assumptions, such as 
assuming proliferation only occurs near the aggregate ends) could 
easily be incorporated. 

Several authors [8,10,24,25] have used Monod kinetics to 
describe the relationship between oxygen concentration and cell 



proliferation, so that the proliferation rate is approximately 
constant at low concentrations, increases with concentration for 
intermediate concentrations, and plateaus to a maximum at high 
concentrations. Following McElwain and Ponzo [26], we approx- 
imate Monod kinetics in a piecewise linear fashion as shown in 
Figure 3A, and assume that cells die below a certain oxygen 
concentration. Following O'Dea et al. [7], we take the cell 
proliferation/death rate to be a stepped function of the fluid shear 
stress on the cells along the membrane outer surface. 



\y=H (see Section A. 3 in File SI for 



details on how this is approximated in the reduced model), with 
faster cell proliferation at intermediate shear stresses and cell death 
at high shear stresses as shown in Figure 3B. This relationship is 
based on experimental observations (see Table S6 and references 
therein). 

If there are N cell aggregates at time t, such that the yth 
aggregate occupies the interval [x2/_i(0>^2/(0] (/ = 1^ ■ ■ • then 
the evolution equation for the length of the Jth aggregate is 



dt 



(X2J-X2j-l)-- 



■.v-2/0 



(Gp{Ce(x,Hm, t),(Te^y{x,H„,t))) 



- Gd(ji:e{x,Hm,t),Ce,xy{x,H,„,t)))dx, (22) 

where Gp and Gd are the cell proliferation rate and death rate 
respectively at each point in the aggregate due to the combined 
effects of the oxygen concentration and shear stress (see 
Figure 3C). The cell proliferation and death functions are 
integrated over the length of the aggregate (from X2/-1 to X2j) to 
give the net rate of change of its length. 

There is littie quantitative data available from which to pose 
and parameterise a growth law, so we use the following 
constitutive relationships for the proliferation and death rates to 
capture the key known features of how oxygen and shear stress 
impact growth 

Gp{Ce,(ye,xy) = ((Q - Ci)H(q - Ci)H(C2 - 

+ (C2 - Ci)H(q - C2)) X (^iH(2i - a.^xy) (22a) 

-I- A2li{<7e^y - Sl)H(i;2 - (Te,xy)), 



GdiCe>Ce,xy) = 

Bo{(l+ Bi H(Co - Ce))(l + B2Hiae^y - S2)) - 1) . 



(22b) 



In (22a) and (22b), H is the Heaviside function; Co the oxygen 
concentration below which cells die; Ci the minimum concentra- 
tion for cell proliferation above which the aggregate growth rate 
increases linearly with Ce\y=H„! C2 the concentration above which 
the aggregate growth rate is a constant maximum; Si the shear 
stress at which the transition in the aggregate growth rate from the 
constant baseline rate Ai to the elevated rate A2 occurs; S2 the 
shear stress above which cells die; and Bo, B\, B2 constants that 
together determine the constant death (aggregate shortening) rates 
for low concentration {Cely^jj^ < Cq), high shear stress 
{^e,xy\y = H„>'^2) and a combination of the two (Ce\y=H„<Co, 
ae,xy\y=H„ >£2). These forms of Gp and Gj capture the fact that 
most cells die if they receive insufficient oxygen, remain quiescent 
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ABC Aggregate growth rate (cm/day) 




Figure 3. Prescribed variation in aggregate growth rate with oxygen concentration and shear stress. A, Variation with oxygen 
concentration, cdy^jj^, for fixed shear stress, (t,,,j||.^h_^ =0.15 Pa. B, Variation with shear stress, (t,,wI,=h„,. for fixed oxygen concentration, 
Ce\y=H„ =0.19molm~^. C, Variation with both oxygen concentration and shear stress. Concentration thresholds: Co = 0.01 molm^', C\ =0.05niolni^' 
and C2 = 0.18molm^^. Shear stress thresholds: Ei=0.1Pa, E2 = 0.2Pa. Growth rate parameters: =0.68{molm"')^'day^', A2 = 
1.35 (mol in--'')"' day-', 5o =0.075 day-', 5i =^2 = 1- 
doi:1 0.1 371 /journal.pone.01 0581 3.g003 



in low oxygen conditions, proliferate at an increasing rate with 
increasing oxygen concentration and die if the fluid shear stress 
they experience is too high. Despite the lack of experimental data 
for parameterisation, the model can be used to test the impact of 
promoting either sensitivity to oxygen or shear stress, through the 
concentration and shear stress thresholds and growth and death 
rate constants. 

We make the further assumption that cell division forces cells 
outwards in such a way that the growth at the ends of the 
aggregate is equal, i.e. 

dt dt 2J^^/_i(0 (23) 

- Gd{Ce(,xMn„t),ae,xyi.X,Hn„t))) dx, 

except when one end of the aggregate reaches the end of the 
membrane, in which case there is only growth at the opposite end 
and the total growth is halved. This gives 2N equations for the 2N 
unknown positions of the aggregate ends (/= 1, . . . ,2N). To 
close equations (l)-(23) we prescribe the initial distribution of 
aggregates via x,(0) (/=!,... ,2N). 

Reduced Model 

Typical values of the model parameters, where possible taken 
from the literature or determined by testing in the laboratory, are 
given in Table 1. Equations (l)-(23) are nondimensionalised and 
the small aspect ratio of the fibre lumen [e = H / L = 2 x 10--^« 1), 
membrane and ECS exploited to simplify the system (see Section 
A in FUe SI for details). The bounds of applicability of the reduced 
model are set by the values of the key dimensionless parameters — 
the reduced Reynolds number, reduced Peclet numbers and 
second Damkohler number. The reduced Reynolds number 
measures the ratio of inertial to viscous forces in the fluid 
accounting for the small aspect ratio of the lumen, and is given by 
e^Ke = e^pUL/ /.i, where p is the fluid density and U is the typical 
lumen fluid velocity set by the inlet flow rate 2,„. For typical inlet 
flow rates used in experiments, 0.02-2 ml min-' [1,19], the 
reduced Reynolds number e^Re lies in the range 6x10-''- 
5 X IQ-^ (with p= 1000 kgm-', the density of water at 20°C to 
2 s.f [27]), which justifies our neglect of inertia in the fluid flow 
equations. The reduced Peclet numbers in the different regions. 



e^PCi = e^ UL/Di (i = l,m,e), represent the ratio of the rate of axial 
advection of oxygen (in the x-direction) to radial diffusion (in the 
j-direction). For tiie range of flow rates given above, the reduced 
Peclet numbers are between 0.2 and 177 (see Table SI), so we 
retain terms at order e^Pe, in the leading order oxygen transport 
equations. Finally, the second Damkohler number, 
y = HT/{DeCi„), gives the ratio of the rate of oxygen uptake by 
the cells to the rate of diffusive transport. The jump in oxygen flux 
across the aggregates is assumed to balance the oxygen uptake by 
the cells so y must be of order 1 for the reduced model to be valid. 

Many other model parameters vary with cell type: the inlet 
oxygen concentration, Ci„; the half maximal uptake flux concen- 
tration, C1/2; the concentration and shear stress thresholds for the 
aggregate growth, Co, Ci, C2, 2i, £2; and the growth rate 
parameters, Ai, A2, B\, B2, Bt,. Experimental data for maximal 
oxygen uptake rates and Ci,,, C1/2 and Co values for different cell 
types is presented in Tables S2 and S4 respectively. In the 
aggregate growth simulations in the Results section we foUow 
Shipley et al. [9] and use values of the transport and uptake 
parameters appropriate to rat cardiomyocytes from the study of 
Radisic et al. [28] (see Table 1). We convert Radisic et al.'s 
estimate of the maximal volumetric oxygen uptake rate of the rat 
cardiomyocytes to a value for the maximal oxygen uptake flux, F, 
as follows 

r = uptake rate per cell x surface density of cells on membrane 
= maximal volumetric uptake flux x cell volume 

X surface density of cells on membrane 
= 2.64xl0-^molni-^s-' x LexlO^'m' x 7.3x lO'^cellsm^^ 
= 3.08 X 10" Vol m^^s"', 

where the surface density of cells on the membrane is taken from 
[29] and the average cell volume estimated from the data in Table 
S3 (see Table S2 for further details). 

Since experimental data for the concentration and shear stress 
thresholds is sparse (see Tables S5 and S6 for data for different 
shear stress regimes), a range of values is used in the simulations to 
determine their effect on the aggregate growth. By comparing 
predictions for aggregate growth from these simulations with 
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Table 1. Fluid and oxygen transport and aggregate growth parameter values. 





Parameter 


Description 


Typical Value(s) 


Reference 


H 


lumen radius 


200 iim 


[19] 


H„ 


fibre outer radius 


200 jim 


[19] 


Hu 


ECS outer radius 


500 


[19] 


L 


fibre length 


10 cm 


[19] 


M 


dynamic viscosity of fluid 


10"^Pas 


A 
V 


^ 


membrane porosity 


0.77 


[31] 




membrane permeability 


2.39x lO-^^m^ 


* 


Q,n 


2D lumen inlet flow rate 


3x 10-^-3x lO-^m^s"' 


[19] 


U=QinlH 


typical flow velocity 


1.5x 10^-^-0. 13ms-* 




Pl.au, 


lumen outlet pressure 


1.03-1.31 X 10^ Pa 


★ * 


Kh 


permeability of membrane outer 
surface + cell aggregate 


8x 10-^(kgm-2s-*)-^ 




Km 


membrane outer surface permeability 


8x 10-^(kgm-2s-*)"^ 




D, 


lumen oxygen diffusivity 


3x lO-^m^s-i 


[32] 


D„ 


membrane oxygen diffusivity 


3xl0-'0m2s-' 


[32] 


Du 


ECS oxygen diffusivity 


3x lO-^m^s-* 


[32] 


r 


maximal oxygen uptake flux 


3.08 X lO-'^molm-^s"* 


[28] t 


Cm 


inlet oxygen concentration 


0.22 molm"^ 


[28] t 


Cl/2 


half maximal uptake rate concentration 


6.9x 10--^molm-^ 


[28] t 


Co 


minimum oxygen concentration 
for cell survival 


6x lO-^molm"^ 


[28] t 


Ci 


minimum oxygen concentration for 
cell proliferation 


% 


t 


Ci 


maximal growth rate oxygen concentration 


% 


% 


Si 


shear stress threshold for elevated 
proliferation rate 


+ 


t 


S2 


shear stress threshold for cell death 


t 


t 


J. 


baseline growth rate for o^, ^^Y\y^ ^^^^ "^^i 


u.68{molm ) day 


m 


A2 


growth rate for Zi <(Tt-,A-v v=//,„ 


1.35(molm-')"'day-' 




Bo 


baseline aggregate shortening rate 

for cvlv=//,„ <*^o or cr^-av v=//,„ >^2 


0.075 day-' 




Bx 


weight factor for aggregate 
shortening rate for Ce\y^H„ <Cq 


1 




Bi 


weight factor for aggregate 
shortening rate for fj,, vi lj, ^ >S2 


1 





^ Value for water at 20"C {from http://physics.info/viscosity/}. 

* In all simulations /*/.„„, is chosen to ensure that there is no backflow at the lumen outlet. 

* Determined in experiments performed by L.A.C. Chapman, R.J. Shipley and M.J. Ellis at the Chemical Engineering Department, University of Bath, 
t Values stated are for rat cardiomyocytes (see Tables S2 and S4 for values for other cell types}. 

% Parameter values for which there is limited experimental evidence (see Table S6 for data on shear stress thresholds} and that are varied in the aggregate growth 
simulations. 

doi:l 0.1 371/journal.pone.Cl 0581 S.tOOl 



experimental results, the relative oxygen concentration and sliear 
stress sensitivity of tlie cells can be inferred. Alternatively, 
experiments can be performed to measure the thresholds for a 
particular cell type and the model used to predict how the 
aggregates will grow. The baseline aggregate growth rate, ^i, is 
estimated from measurements of population doubling times of 
human bone marrow stromal cell culture in a HFB [2] (see Section 
B.2.2 in File SI for details). 



Numerical Methods 

The reduced system of equations for the flow, oxygen transport 
and cell aggregate growth is solved using code developed in the 
numerical programming software MATLAB (distributed by 
Mathworks Inc., see http://www.mathworks.co.uk for details). A 
description of the algorithm is given in Section C.3 of File SI. The 
code was validated by comparing numerical solutions of the 
reduced model without aggregate growth to solutions of the fuU 
fluid and oxygen transport model (given by equations (1)-(16)) 
computed using the finite element software COMSOL Multi- 
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Figure 4. Aggregate distribution witKi oxygen and shear stress distributions at different time points in a growtli simulation. A-C, 

Aggregate distribution with oxygen distribution. D-F, Aggregate distribution with shear stress distribution. Aggregates shown in white. Solid black 
lines are flow streamlines. The dark blue in the lumen and membrane in D-F indicates very low and negative fluid shear stress in these regions 
respectively. Parameter values: ^ = 5, L„j,,j, = 0.5cm, e,„ = 3 x IQ-^m^s"', P;,<,„, = 1.31 x lO'Pa, r = 3.08 x IQ-'molm-^s-', Ci/2 = 
6.9x lO-^molm-3, Co = 6x 10--^molm-^ Q =0.05molm-\ C2 = 0.18molm-3, Si =6.03 Pa, = 0.08 Pa, A2 = 2Ai, BQ=AtQ„/2, Bi=B2 = l. 
All other parameter values as in Table 1. 
doi:1 0.1 371 /journal.pone.01 0581 3.g004 



physics 4.3 (distributed by COMSOL Inc., see http://www. 
comsoI.com for details). Over the range of parameter values for 
which the approximation in the reduced model is valid, the error 
in the reduced model solutions was consistent with the size of 
terms neglected (results not shown). 



Results 

Example Oxygen, Shear Stress and Aggregate 
Distributions 

Figure 4 shows the distribution of rat cardiomyocyte aggregates 
along the membrane with the corresponding oxygen and shear 
stress distributions in the ECS at different time points in a growth 




X (m) X (m) 

Figure 5. Advective oxygen flux across membrane outer surface for different flow rates and outlet pressures for rat cardiomyocyte 
aggregates. A, Advective flux, Ceiv|,,^H , for different flow rates. Arrow shows direction of increasing Qj„ (Pi,out = 1.03 x 10'' Pa). B, Flux for different 
outlet pressures. Arrow shows direction of increasing Pi^om (Qm = ^ x 10^'m^s^'). Initial aggregate distribution as in Figure 4A. Parameter values: 
Ci =0.1molm"^ C2 = 0.21 molm"\ Ei = 1 Pa, E2 = 2Pa, A2 = 0, Bo=AiCi„/2, Bi = 1, ^2 = 0. All other parameter values as in Figure 4. 
doi:1 0.1 371/journal.pone.01 0581 3.g005 
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Figure 6. Shear stress along membrane for different flow rates and outlet presures. A, Shear stress for different flow rates. Arrow shows 
direction of increasing g;„ (/"/u,,, = 1.03 x lO^Pa). B, Shear stress for different outlet pressures. Arrow shows direction of increasing P/,„„, 
(£?m = 3 X All other parameter values as in Figure 5. 

doi:1 0.1 371 /journal.pone.01 0581 3.g006 



simulation for candidate parameter values (see figure legend). The 
initial aggregate distribution (shown in Figures 4A,D) consists of 
five aggregates of length 0.5 cm distributed evenly along the 
membrane. It is used for all subsequent simulations, unless 
otherwise stated, so that the impact of oxygen versus fluid shear 
stress on growth can be isolated and compared. Due to uptake by 
the cells, oxygen is depleted in regions that surround each 
aggregate and extend outwards from the membrane into the ECS. 
These regions increase in size as the aggregates grow and 
eventually merge to form one large region (Figures 4A-C). The 
flow of the culture medium ensures that the oxygen concentration 
remains highest at the upstream end of the ECS until the 
aggregates come into contact with it. The shear stress profile along 
the membrane remains fairly constant in time as the aggregates 
grow (Figures 4D-F), because the aggregates have only a small 
effect on the membrane outer surface permeability (for the chosen 
values of Kio and Khi). The aggregates do not reach confluence 
after 40 days (the longest culture period reported in the literature 
[30]) as the shear stress exceeds the cell death threshold 
{'^e,xy\y=H,„>^2) towards the outlet. Thus the three aggregates 
nearest the inlet grow and coalesce to form one large aggregate, 
while the fourth aggregate remains static and the rightmost one 
shrinks. 

Influence of Flow on Aggregate Growth 

Impact of Flow on Oxygen Flux and Shear Stress. To 

understand how the inlet flow rate, Qj„, and outlet pressure, P/,o,«, 
affect the growth of the aggregates it is necessary to determine how 
they affect the oxygen concentration around, and shear stress on, 
the aggregates. 

Figures 5A and 5B show how the advective oxygen flux through 
the outer surface of the membrane, QV^j^^^^, changes when 2,„ 
and P/,oi// are varied for the static rat cardiomyocyte aggregate 
distribution. The flux dips in the regions of the membrane covered 
by aggregates because of the oxygen uptake and reduced fluid 
permeability there. Increasing g,„ or Pi^oui increases the oxygen 
flux through the membrane. However, the system is more sensitive 
to changes in Pi^oui than Qi„: a 100-fold increase in Qi„ (from 
3xl0^'m^s^' to 3xI0^'m^s^') leads to an increase of at 
most 50% in the flux, while a 14-fold increase in (Pi^mii~ Patm) 



(the outiet pressure relative to atmospheric pressure, from 
2.I4kPa to 29.7kPa, or 0.31 to 4.31 psi) produces an 8-fold 
increase. This is because the bulk of the pressure drop from the 
lumen inlet to the ECS outlet occurs across the membrane, so 
increasing Pi^out gives a much bigger pressure gradient across, and 
therefore flow through, the membrane. Similarly, the shear stress, 
ae^xy\y=H,„^ is found to increase with distance along the membrane 
(from 0 at the upstream end, which is closed) and as Qin and Pi,out 
are increased, the effect being more pronounced for changes in 
Pi, out than Qin (see Figure 6). The 100-fold increase in Qi„ 
produces a 26% increase in the maximum shear stress over the 
membrane, whereas the 14-fold increase in (Pi,out — Paim) causes a 
10-foId increase in the maximum shear stress. 

Impact of Flow on Aggregate Growth. First we distinguish 
the effects of oxygen and shear stress on aggregate growth as Qi„ 
and Pj^out are varied by comparing simulations for which growth is 
purely oxygen-dependent with those for which it is purely shear- 
stress-dependent. Oxygen-dependent growth is appropriate to cell 
types such as endothelial cells and hepatocytes, which are tolerant 
to high shear stresses or have high oxygen demands, while shear- 
stress-dependent growth is appropriate to cell types such as 
chondrocytes and mouse embryonic stem cells, which have low 
oxygen demands or are sensitive to changes in shear stress (see 
Tables S2 and S6). In the model, purely oxygen-dependent growth 
is achieved by setting the shear stress growth parameters 
^2 = ^2 = 0 in Equation (23) and choosing S] such that 
ae,xy\y=H„, always; and purely shear-stress-dependent growth 
is achieved by taking B2 = Q m Equation (23) and C2 low enough 
such that Ce\y^jj„,>C2 always. Following this, we consider 
aggregate growth for cells that are sensitive to both oxygen levels 
and shear stress, choosing the concentration thresholds 
C\ =0.05 molm^' and C2 =0.18 molni^' in Equation (23) such 
that the growth rate is linearly proportional to the concentration 
over a wide range and using low values for the shear stress 
thresholds 2i =0.01 Pa and £2 = 0.05 Pa. Oxygen- and shear- 
stress-dependent growth is relevant to cell types such as 
cardiomyocytes, mesenchymal stem cells and fibroblasts (see 
Tables S2 and S6). 

Simulations are run until either the aggregates reach confluence 
or 40 days have elapsed. In the former case the confluence time, tc. 
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Figure 7. Relationship between confluence time and outlet pressure for different flow rates for rat cardiomyocyte aggregates. A, 

Relationship for oxygen-concentration-dependent growth. B, Relationship for shear-stress-dependent growth. Time step: Ai = 3.24hrs 
(dimensionless time step 10^-). Parameter values in all as in Figure 5. Parameter values in B: Ci =0.02molm 

C2 = 0.08molm"^ Ii=0.05Pa, E2 = 0.16Pa, ^2 = 2^1, Bo=A,Ci„/2, Bi=0, £2 = 2. All other parameter values as in Figure 4. Curves stop as 
backflow occurs at higher outlet pressures. 
doi:1 0.1 371 /journal.pone.01 0581 3.g007 



is recorded, in the latter case tlie total lengtli of the aggregates, 
Lioi, is recorded (at confluence L,oi = L= 10cm). 

Oxygen-dependent Growth. Figure 7A shows how tlie 
confluence time, /f, varies with Q/„ and Pi^out for simulated rat 
cardiomyocyte aggregate growth. As Pi^out increases with Q„, 
fixed, the curves end at the maximum oudet pressure for which 
there is no backflow. As expected, ?c decreases as Qifj and Pi^oui 
increase, the effect being more marked for changes in Pi^oui- This is 
because increasing Qi„ or P/_o„, increases the advective flux of 
oxygen through the membrane and therefore the oxygen 
concentration around the aggregates and their growth rate. 

Shear-stress-dependent Growth. The results presented in 
Figure 7B reveal that as Pi^out increases there is a sharp decrease in 



the confluence time at around -P/,ou/ = 1-09 x 10"" Pa (15.8psi), but 
that changes in Qj„ have litde effect on tc- As Pi^out increases 
through Pioui = 1 -09 x lO' Pa, it is observed in the simulations that 
the number of aggregates on which the shear stress is higher than 
S] increases sharply, leading to a significant increase in the net 
aggregate growth. The effect is barely noticeable when 2„, is 
increased, because the associated increase in the fluid flux through 
the membrane is much smaller. 

Combined Oxygen-concentration- and Shear-stress- 
dependent Growth. Figure 8A shows that aggregates of rat 
cardiomyocytes with the oxygen and shear stress sensitivities given 
above only reach confluence when Pj^out 

< 1.09 X 10^ Pa, and that 
for Pi^oui> 1-09 X 10^ Pa, L,oi decreases as Pi^oui increases, up to 
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Figure 8. Oxygen-concentration- and shear-stress-dependent growth of shear-sensitive rat cardiomyocytes. A, Relationship between 
final total aggregate length and outlet pressure for different flow rates (graphs overlap as there is no change with flow rate). B, Variation in 
confluence time with flow rate and outlet pressure for < 1.09 x 10' Pa. Parameter values: Q =0.05molm^', C2 = 0.18 molm^^, Zi =0.01 Pa, 
E2 = 0.05 Pa, A2 = 2Ai, Bo=A\Cj„l2, B\=B2 = 1. All other parameter values as in Figure 4. 
doi:1 0.1 371/journai.pone.01 0581 3.g008 
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Figure 9. Aggregate growth and oxygen distribution at different time points for different initial aggregate distributions and low 
outlet pressure. First row: left-skewed distribution. Second row: right-skewed distribution. Third row: even distribution. Solid black lines are flow 
streamlines. Parameter values: N = 5, L„gj; = 0.5cm, Q,„ = 3 x IQ-'m^s-', Cj =0.05molm-^ C2 = 0.18molm-^ Ei=0.03Pa, i:2=0.08Pa, 
A2 = 2A\, B(} = AiCi„/2, Bi=B2 = l. All other parameter values as in Figure 4. 
doi:1 0.1 371/journal.pone.01 0581 3.g009 



1.17x10' Pa (17 psi) (above which it is constant). This is because 
at outlet pressures above 1.09 x 10' Pa the shear stress exceeds the 
cell death threshold ^2 over more of the membrane and this effect 
outweighs the positive effect of greater oxygen delivery to the cells. 

Figure 8B shows how the confluence time changes with flow 
rate and outlet pressure for flow rates between 3x lO^^m^s^' 
and 3x10^' w? s ^ ' and outlet pressures ranging from 
1.035xl0'Pa to 1.09xlO'Pa (or 15psi to 15.8psi) (ranges 
achievable with standard laboratory equipment). As P/^oui 
increases, the confluence time first decreases, for P/^our up to 
1.083 X 10' Pa (15.7psi), and then increases. This is because the 
region of the membrane in which the proliferation rate is elevated 
(i.e. where 2i <(T(.j; |v=i/,„ <2^2) increases in extent, covering 
more of the aggregates, as Pi^oui increases up to 
■P/,oi/r = 1-083 X 10' Pa, but the region in which the shear stress 
exceeds the death threshold £2 becomes larger for 
P/,o„,>1.083x lO'Pa. For P/,,„„ < 1.05 x lO'Pa (15.2psi), the 



confluence time decreases marginally as Qi„ increases due to 
greater oxygen flux through the membrane. In the range 
1.05 x 10'Pa<P/,o„,< 1.07 x lO'Pa (15.2psi<P/,o„,<15.5psi), 
the conffuence time remains approximately constant as Qi„ 
increases. For 1.07 x lO' Pa <P/,o„,< 1.083 x lO' Pa, the conffu- 
ence time decreases as Qi„ increases, reaching a minimum of 7.9 
days for 2.75 x IQ-'m^ s"' < g,,, < 3 x lO-'m^s-' and 1.077 x 
10'Pa<P/ o„,<1.083x lO'Pa. For 1.083 x 10'Pa<P/,<,„,< 
1.09 X 10' Pa the confluence time increases with increasing Qi„, 
since the detrimental effects of the shear stress exceeding S2 start 
to outweigh those of increased oxygen delivery. Hence, 
2.75x 10-'m2s-i<e,„<3x lO-'m^s-i and 1.077 x 10' Pa< 
P/,oM/ < 1 ■ 08 3 X 1 0' Pa ( 1 5 . 6 psi < P/,o„, < 1 5 . 7 psi) are predicted to 
be the optimal ffow rate and oudet pressure ranges for 
proliferation of these shear-sensitive rat cardiomyocytes. 
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Influence of Cell Seeding on Growth 

Although the model can be used to simulate the growth of real 
seeding distributions, since seeding data is not available we use 
three different simple initial distributions to test the impact of skew 
in the distribution (towards the inlet end or outlet end) on 
confluence time. We expect the skew to have an effect because the 
oxygen and shear stress distributions that arise are spatially non- 
uniform and evolve as the aggregates grow. Each of the 
distributions has five aggregates each of initial length 0.5 cm, 
and the aggregate growth is simulated for a low outlet pressure 
(-P/,oM/ = 1-03 X lO^Pa) and a high outlet pressure (Pi,oui = 
1.31xlO'Pa). (We do not vary Qi„ since its effect on the 
transmembrane fluid flux is much weaker.) The initial distributions 
we consider are: (i) left-skewed (i.e. aggregates concentrated 
towards the inlet); (ii) right-skewed (i.e. aggregates concentrated 
towards the outlet); and (iii) even (i.e. aggregates spaced uniformly 
across the membrane) (see / = 0 panels in Figure 9). 

With the lower outlet pressure the aggregates grow to 
confluence in all three cases, albeit at differing rates. The 
aggregates reach confluence fastest if they are right-skewed 
(?£ = 10.5 days), take longer if they are evenly distributed 
(/c = 17 days), and take more than twice as long (?(. = 24.8 days) 
if they are left-skewed. For the right-skewed distribution the 
aggregates start further away from, and grow towards, the inlet. 
Hence, once the aggregates have merged to form one large 
aggregate, the end of the aggregate nearest the inlet continues to 
experience a high oxygen concentration as it grows, despite the 



weak oxygen flux through the membrane from the low oudet 
pressure (Figure 9). Changes in the shear stress along tiie 
membrane with the initial aggregate distribution do not have 
any impact on the rate of growth to confluence as the shear stress 
is always below the threshold for increased proliferation 
(Si =0.03 Pa) for all four distributions (results not shown). 

With the higher outlet pressure, the aggregates fail to reach 
confluence after 40 days for all four initial distributions. The total 
aggregate length at 40 days is greatest for the even distribution 
(7 cm) and slightiy less for the left-skewed and right-skewed 
distributions (5.76 cm and 5.84cm respectively). This is because 
when the outiet pressure is large the oxygen concentration around 
the aggregates is similar for the difiFerent distributions but the shear 
stress exceeds the threshold for cell death near the outlet end of the 
membrane. Aggregates at the oudet end therefore tend to shrink, 
so initial differences in growth between the distributions disappear 
over time and the aggregates merge to form one large aggregate at 
the inlet. 

Discussion 

We have considered how the lumen inlet flow rate iQ,„ and 
outlet pressure Pi^out affect the growth of cell aggregates with 
different oxygen and shear stress sensitivites through their impact 
on the oxygen delivery to and shear stress on the cells. Increasing 
Qin or Pi^out increases the transmembrane fluid flux, the increase 
being more pronounced for increments in Pi^om over the range of 
flow rates and pressures typically used in experiments. The 



A 




T- 10 '' - 0.05 Pa 0.05 - 0.2 Pa > 0.2 Pa 



Increasing shear stress tolerance 

Figure 10. Optimal HFB operating conditions for proliferation of cell types witKi different oxygen uptake rates and shear stress 
tolerances. Wedge width and shading show relative values of flow rate, 2„„ and outlet pressure, /"/_„„/, that give maximum cell yield for cells with 
different oxygen demands and shear tolerance; greater widths and darker shades representing greater values, up to g„, = 3 x 10^'m-s^' and 
Pi.oiii = 1.31 X 10'' Pa. Angle of slanted line in crossed lines shows direction and degree of skew in position of initial aggregate distribution that gives 
maximum cell yield with optimal and g„, values; right slanting representing skew towards outlet and left slanting skew towards inlet end and a 
greater angle representing greater skew. Box marks optimal operating conditions for cells with oxygen uptake rate of ~ 10^'^mols^' cell^' and 
shear tolerance of E2 ~0.1 Pa (see main text). 
doi:1 0.1 371/journal.pone.01 0581 S.gOlO 
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increase in transmembrane fluid flux tias two eflects: more oxygen 
is adverted through the membrane to the aggregates, and the 
horizontal flow in the ECS is stronger so the aggregates experience 
higher shear stress. The higher oxygen concentration increases the 
rate of aggregate growth, but the eflfect of the increased shear stress 
depends on the shear-sensitivity of the cells (i.e. the shear 
thresholds 2i and S2, which mark the transitions to elevated 
proliferation rate and cell death) and the aggregate distribution, 
since the shear stress is non-uniform over the membrane. The 
shear stress increases from zero at the inlet end of the membrane 
to a maximum at the outlet end, which ranges between 0.01 Pa 
and 0.2 Pa for the range of flow rates and outlet pressures 
considered (3 x IQ-'^m^s"' < 2,„ <3 x IQ-^m^s"' and 1.03 x 
10^Pa<P/,<,„,<1.31 X lO^Pa (15psi<P/,„„,< 19psi)). Figure 10 
summarises the operating conditions and basic seeding strategy 
that maximise the cell yield for cells with different oxygen uptake 
rates and shear stress tolerances. The optimal initial position of the 
aggregates is described in terms of the skew of the aggregate 
distribution towards the ECS outiet. For rat cardiomyocytes with 
an uptake rate of 4.22x 10^'^mol s^' ceU^' and shear tolerance 
of S2 = 0.05 Pa, for example, the optimal operating conditions and 
aggregate seeding strategy are 2.75x 10-5m2s-'<e,„< 
3xl0-5m2s-' and 1.077 x 105Pa<P/,„„,< 1.083 x lO^Pa 
(1 5.6 psi < P/ o„, < 15.7 psi) and an even distribution. 

It is evident from Figure 10 that achieving the maximum cell 
yi(;ld recjuires a careful balance of the flow rate, outlet pressure and 
seeding distribution. For instance, the optimal seeding distribution 
of the aggregates changes from moving towards the outlet with 
increasing oxygen demand when S2 > 0.2 Pa to moving towards 
the inlet when S2 drops below 0.05 Pa. This is because the low 
slu-ar tolcranix- of tlu- cells prevents aggregates from growing near 
the outlet despite the improved oxygen delivery that can be 
achieved by seeding them there. The extent to which the ceU 
seeding process can be controlled is thus an important issue. 
Indeed, the more random it is, the more likely aggregates are to 
form in regions with unfavorable growth conditions. This 
motivates further experimental studies to improve the control 
over the initial cell distribution achievable with seeding protocols. 
In the meantime, for random seeding of cells with typical oxygen 
demands, the cell yield can be optimised by using a flow rate and 
outiet pressure that ensure the shear stress is in the elevated 
proliferation rate range over as much of the membrane as possible. 
For example, from Figure 10 (se(^ box) wc can see that for cells 
with an oxygen uptake rate of the order of 10^'^mols^' cell^' 
and a shear stress death threshold 22~0.1Pa in an initial 
aggregate distribution skewed slightly towards the outlet, a flow 
rate of approximately 3 x lO^^m^s^' and an outlet pressure of 
1.10 x 10^ Pa (16psi) will ensure that the region of shear-stress- 
enhanced proliferation is as large as possible. 
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